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Autonomous  Distant  Visual  Surveillance  of  Satellites 

John  E.  Mclnroy,  Senior  Member,  IEEE,  Lawrence  M.  Robertson  and  R.  Scott  Erwin 


Abstract —  This  paper  develops  three  new,  interconnected  tech¬ 
niques  useful  for  the  autonomous  distant  visual  inspection  of 
satellites.  First,  silhouetting  of  man  made,  erratically  illuminated 
satellites  is  performed.  Illumination  cases  include  full  sun  from 
an  arbitrary  (often  awkward)  viewing  angle  and  un-illuminated 
(back-lit  by  the  star  field).  New  Statistical  Straight  Line  Snakes 
(SSLS)  prove  efficient  in  finding  the  silhouette,  even  in  the  un¬ 
illuminated  case.  The  silhouette  is  composed  of  straight  line 
segments,  which  are  easy  to  calculate,  fit  the  straight  lines 
inherent  in  man  made  objects,  and  lend  themselves  to  further 
processing  (pose  estimation,  template  matching,  etc.)  Once  the 
silhouette  has  been  used  to  find  correspondence  points,  a  second 
method  for  detecting  a  moving,  nearby  chaser  vehicle  is  derived. 
The  hard  case  is  treated  in  which  the  chaser  and  satellite  are  so 
nearby  that  their  images  are  blurred  together.  The  algorithm 
finds  the  dimension  of  motion  generated  by  the  sequence  of 
images.  If  the  dimension  is  higher  than  that  explained  by  a  single 
rigid  body,  then  this  indicates  a  possible  chaser.  No  knowledge  of 
the  satellite  or  chaser’s  shape  or  motion  is  required.  Independent 
relative  motion  between  the  satellite  and  chaser  is  required- 
if  the  chaser  is  immobile  with  respect  to  the  satellite,  then  a 
third  technique  must  be  used.  This  third  method  incorporates  the 
satellite’s  solid  model  to  estimate  its  pose  from  a  noisy,  diffraction 
limited  image.  The  pose  is  then  combined  with  the  solid  and 
optical  model  to  create  synthetic  expected  images.  Inspection 
is  performed  by  comparing  these  with  the  actual  images.  The 
new  pose  algorithm  first  estimates  depth  by  a  least  upper  bound 
technique.  A  fast  method  is  derived  of  optimally  estimating  the 
rotation  matrix  by  a  sequence  of  analytical  solutions  (rather  than 
a  nonlinear  numerical  optimization!).  Simulations  illustrate  the 
use  of  all  three  techniques  on  images  obtained  when  viewing  low 
Earth  orbit  satellites  from  the  ground. 

Index  Terms —  affine  cameras,  pose  estimation,  visual  servoing, 
template  matching,  object  recognition,  motion  subspace 


I.  Introduction 

Visual  surveillance  of  satellites  presents  new  challenges  and 
opportunities  for  computer  vision.  The  distance  between  an 
Earth  based  telescope  and  a  satellite  is  great,  which  implies 
that  the  images  are  often  blurry  due  to  the  effects  of  the 
atmosphere,  diffraction,  etc.  Moreover,  illumination  is  erratic 
and  harsh,  varying  from  a  very  sharp  contrast,  to  partial 
illumination,  to  un-illuminated  but  back- lit  by  stars.  Finally, 
the  objectives  are  somewhat  different.  For  instance,  it  may 
be  of  interest  just  to  detect  that  any  body  is  in  the  image. 
Alternatively,  it  may  be  desired  to  know  if  a  single  blurry 
image  contains  multiple  bodies. 

This  work  was  supported  by  the  National  Research  Council  and  Air  Force 
Office  of  Scientific  Research. 
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This  paper  develops  three  new,  complementary  tools  for  the 
distant  visual  surveillance  of  satellites.  First,  silhouetting  of 
man  made,  erratically  illuminated  space  objects  will  be  per¬ 
formed.  Second,  methods  for  detecting  multiple  independently 
moving,  visually  overlaid  bodies  will  be  found.  Third,  methods 
for  inspecting  known  satellites  will  be  derived. 

Finding  the  silhouette  of  spaceborne  objects  is  useful  for 
both  detecting  the  presence  of  an  object  and  for  finding 
correspondences  between  images  and  models.  When  viewing 
from  Earth,  sun  illumination  plays  a  significant  role,  as  it 
can  be  very  harsh,  partial,  or  even  completely  missing.  The 
un-illuminated  case  is  particularly  unique,  as  much  can  still 
be  done  due  to  back  lighting  by  stars.  Traditional  gradient 
based  edge  detectors  ([9],  [10])  fail  utterly  in  the  unilluminated 
case.  More  recent  contour  growing  (’’snake”)  techniques  ([17], 
[14],  [1],  [3])  are  more  promising.  This  paper  develops  a 
new  Statistical  Straight  Line  Snake  (SSLS)  algorithm  which 
is  simple  and  naturally  fits  man  made  objects  with  their 
preponderance  of  straight  edges.  It  directly  outputs  lines  and 
vertices  ready  for  use  in  correspondence,  pose,  and  template 
matching  algorithms. 

Nearby  moving  objects  are  always  a  concern  for  satellites. 
Unfortunately,  visual  evaluation  from  Earth  can  be  difficult 
due  to  blurry  images  caused  by  diffraction,  atmospheric, 
and  scaling  effects.  Satellites  are  in  constant  motion:  is  the 
motion  seen  in  an  image  caused  by  a  single  satellite,  or  are 
multiple  bodies  present?  This  paper  modifies  a  robotic  motion 
segmentation  algorithm  [18]  to  answer  this  question.  No  model 
of  the  motion,  satellite,  or  other  objects  is  required. 

On  the  other  hand,  relative  motion  between  the  multiple 
bodies  is  required.  To  inspect  known  satellites,  and  especially 
to  detect  bodies  immobile  with  respect  to  the  satellite,  a 
new  affine  camera  pose  estimation  algorithm  is  derived.  Once 
the  pose  is  found,  it  is  combined  with  optical  models  and 
solid  satellite  models  to  predict  the  image.  Large  differences 
indicate  an  anomaly. 

Despite  the  availability  of  many  pose  estimation  algorithms 
(see  [11]  [4]  [6]  [15],  [8],  and  [2]  for  starters),  there  remains 
a  need  for  a  rapid,  optimal,  dependable  algorithm  suitable  for 
a  small  number  of  point  correspondences  (but  more  than  n  = 
4).  For  instance,  [16]  recently  added  high  bandwidth  inertial 
sensing  to  complement  the  available  pose  estimation  because 
it  alone  had  too  low  of  a  bandwidth.  For  the  case  of  affine 
imaging,  this  paper  develops  a  new  method  that  is  closed  form, 
yet  provides  optimal  estimates  even  when  n  >  4.  It  does  so 
by  solving  the  standard  3d-3d  optimal  orientation  problem  2n 
times,  where  n  is  the  number  of  data  points.  Since  the  optimal 
orientation  problem  can  be  very  quickly  solved  as  a  4  x  4 
matrix  eigenvalue  problem,  the  overall  computational  burden 
is  light.  A  least  upper  bound  estimate  of  the  scaling  (depth)  is 
proposed  and  combined  with  the  orientation  algorithm  to  find 
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the  complete  pose. 

This  paper  explains  these  techniques  in  roughly  the  same 
order  they  would  typically  be  performed.  First,  a  method  of 
finding  the  silhouette  will  be  described.  It  outputs  point  cor¬ 
respondences  that  are  needed  by  both  of  the  other  algorithms. 
This  data  can  be  directly  used  for  determining  if  the  sequential 
images  contain  multiple  bodies.  Further  inspections  can  be 
performed  by  incorporating  the  satellite’s  known  shape. 

II.  Silhouette  of  Man  Made,  Erratically 
Illuminated  Objects  in  Space 

This  paper  seeks  to  develop  a  technique  for  silhouetting 
which  builds  on  existing  methods,  but  modifies  them  to  be 
especially  robust  when  finding  point  or  line  correspondences 
on  man  made  objects  under  the  erratic  illumination  found  in 
space.  This  may  include  partial  illumination,  zero  illumination 
against  a  star  filled  background,  harsh  illumination,  and  bright 
reflections.  The  method  is  inspired  by  contour  (or  snake) 
growing  ([17],  [14],  [1]),  but  exploits  the  inherent  straight  line 
structure  dominant  in  man  made  objects  such  as  satellites,  as 
well  as  the  specular  nature  of  the  star  filled  background.  Rather 
than  fitting  an  arbitrary  contour  around  an  object,  it  forms  a 
contour  from  line  segments.  The  length  and  number  of  the 
segments  is  iteratively  adapted  to  fit  the  measurements. 

A  brief  outline  of  the  steps  are  as  follows: 

1)  Place  an  initial,  small  polygon  inside  the  object’s  image. 
Store  it  as  a  list  of  vertex  points  proceeding  in  a  counter¬ 
clockwise  direction  around  the  polygon. 

2)  Eliminate  a  vertex  if  the  line  between  adjacent  vertices 
is  shorter  than  the  desired  resolution. 

3)  Add  vertices  to  keep  each  edge  a  length  Id  or  shorter. 

4)  Create  a  set  of  points  along  a  straight  line  edge  between 
vertices.  If  a  point  is  inside  the  object,  then  perturb  that 
point  outward  along  the  edge  normals.  Repeat  this  for 
all  edges  along  the  polygon. 

5)  Fit  lines  to  these  perturbed  points  along  each  edge. 

6)  Create  new  vertices  by  intersecting  non-collinear  adja¬ 
cent  lines. 

7)  Check  for  convergence,  returning  to  step  2  until  conver¬ 
gence  takes  place. 

The  steps  will  now  be  covered  in  more  detail. 

A.  Step  1:  Place  an  initial,  small  polygon  inside  the  object's 
image.  Store  it  as  a  list  of  vertex  points  proceeding  in  a 
counter-clockwise  direction  around  the  polygon. 

An  initial  point  inside  the  object’s  image  can  be  found 
by  creating  a  mask  which  encodes  the  object’s  rough  shape, 
then  moving  it  over  the  image  until  the  best  fit  is  found. 
The  centroid  of  the  mask  giving  the  best  match  is  the  initial 
point.  A  small  circular  polygon  (like  a  hexagon  or  octagon) 
is  then  placed  inside  the  object,  centered  at  the  initial  point. 
The  number  of  sides  in  this  initial  polygon  is  not  especially 
important,  as  vertices  will  be  added  and  removed  during 
the  iterations  to  enforce  edge  length  and  non-collinearity 
specifications. 


B.  Step  2:  Eliminate  a  vertex  if  the  line  between  adjacent 
vertices  is  short. 

Let  dj  and  dj+ 1  (both  in  M2)  denote  two  adjacent  vertices. 
The  length  of  edge  j  can  be  found  from  lj  =  \\dj^g  —  dj\\. 
If  this  length  is  below  the  desired  resolution  (, l min X  then 
vertex  j  +  1  is  dropped  from  the  list  of  vertices.  This  process 
intentionally  limits  the  resolution  of  the  contour  fit,  which  is 
useful,  for  instance,  to  prevent  noise  from  creating  false  edges. 

C.  Step  3:  Add  vertices  to  keep  each  edge  a  length  Id  or 
shorter. 

Let  ^  be  a  desired  target  length  for  each  edge.  If  lj  >  Idi 
then  a  new  vertex  can  be  introduced  between  dj  and  dj+ 1 
at  the  midpoint.  This  is  performed  by  inserting  the  additional 
vertex,  (dj  +  dj+i)/2.  The  length  of  this  segment  (from  dj 
to  the  new  dj+ 1)  can  then  be  evaluated,  and  if  larger  than 
Id,  another  vertex  can  be  added  at  the  midpoint  of  dj  and 
the  newly  created  dj+\.  This  process  is  repeated  until  the 
length  between  dj  and  the  next  vertex  is  less  than  Id,  then 
j  is  incremented.  The  distance  between  the  (possibly  newly 
created)  vertex  j  +  1  and  the  (possibly  newly  created)  vertex 
j  +  2  is  compared  to  Id,  and  so  forth.  The  process  repeats  until 
all  edges  have  a  length  less  than  Id . 

In  the  next  step,  new  edges  will  be  found  by  perturbing 
points  along  each  edge  until  they  reach  the  boundary,  then 
fitting  a  line  to  each  edge  from  these  perturbed  points.  By 
keeping  the  edges  from  growing  too  large,  excessive  averaging 
is  avoided,  thus  allowing  the  contour  to  conform  to  the 
boundary. 

D.  Step  4:  Create  a  set  of  points  along  a  straight  line  edge 
between  vertices.  If  a  point  is  inside  the  object,  then  perturb 
that  point  outward  along  the  edge  normals. 

The  unit  vector  pointing  between  adjacent  vertices  is 

-  _  j-  ~  dj 

I  l^j+i  ~  dj  II 

Points,  x,  along  this  line  segment  can  be  created  by  x  =  dj  + 
olv  where  a  is  between  0  and  the  length  of  the  line  segment,  lj . 
Let  u  denote  the  outward  pointing  normal.  Since  the  vertices 
are  in  a  counter-clockwise  direction,  the  cross  product  u  x  v 
is  positive.  This  implies  that 

r  v2 

u  = 

.  ~V1 

The  points  can  then  be  perturbed  along  the  edge  normal  by 

x  +  (p(I(x)\0)  —  p(I(x)\B))gu  — »  x  (1) 

where  p(I(x)\0)  denotes  the  probability  that  intensity  /(•) 
occurs  at  location  x,  given  that  the  point  lies  on  the  object, 
p(I(x)\B)  denotes  the  probability  that  intensity  /(•)  occurs  at 
location  x,  given  that  the  point  lies  on  the  background,  and  g 
is  a  positive  scalar  gain. 

Since  the  polygon  is  initially  placed  somewhere  deep  inside 
the  object,  initially  it  is  desirable  to  have  a  large  gain,  g , 
that  will  perturb  the  data  points  several  pixels  per  iteration. 
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Once  the  edges  approach  the  boundary,  then  the  gain  should 
be  reduced  so  that  a  fine  convergence  to  the  edges  will  result. 
There  are  many  possible  functions  which  perform  this  action. 
One  method  is  to  make  g  a  function  of  the  iteration  number, 
decreasing  it  linearly  as  the  iterations  proceed.  Another  method 
is  to  make  g  a  function  of  e  =  p(I(x)\0 )  —  p(I(x)\B);  e  is 
between  -1  and  1.  Since  the  polygon  is  initially  inside  the 
object,  e  starts  at  e  =  1.  Once  the  polygon  is  close  to  the 
boundary,  e  — >  —  1.  To  begin  with  a  large  gain,  gmax ,  which 
exponentially  decreases  with  e  and  ends  (when  e  =  —  1)  with 
the  small  gain,  gmin ,  let 

g(e)  =a  +  ebe  (2) 

where 

7  .  1  9  max  9min  h 

b  =  arcsmh - - - ,  a  =  gmax  -  e  (3) 

The  star-filled  backdrop  in  space  images  creates  unique 
problems  as  well  as  opportunities.  When  the  satellite  is  illumi¬ 
nated,  the  main  problem  is  the  stars  act  as  a  source  of  noise, 
making  it  hard  to  tell  (at  distinct  points)  where  the  spacecraft 
ends  and  the  stars  begin.  The  straight  line  averaging  explained 
in  the  next  section  has  proven  robust  to  this  noise. 

On  the  other  hand,  the  satellite  is  back  lit  by  stars,  thus 
it  forms  a  (dark)  image  when  it  is  not  illuminated  by  the 
sun.  This  provides  an  opportunity  to  image  the  silhouette  even 
without  the  sun’s  illumination.  Consider  a  binary  image  where 
pixels  equaling  0  are  black  and  1  are  white.  Under  these 
circumstances,  the  complement  of  the  un-illuminated  satellite 
image  will  be  white.  However,  the  noise  level  will  also  be 
complemented.  If  the  sky  in  the  field  of  view  has  a  fraction 
of  the  total  area  (p  stars)  containing  stars,  then  the  fraction  of 
the  background  containing  noise  in  the  un-illuminated  case  is 
1  —  p stars-  Since  p stars  is  typically  <C  1,  this  creates  a  high 
level  of  noise.  To  handle  this,  a  layer  within  the  boundary  can 
be  checked  and  penalized  for  the  presence  of  noise.  Simply 
check  points  parallel  to  the  edge  between  vertices  j  and  j  + 1, 
but  inside  the  area.  If  a  point  is  believed  to  be  from  the 
background,  rather  than  the  object,  then  perturb  that  point 
inward,  and  incorporate  it  in  the  list  of  line  fit  points.  Let 
pw  be  the  width  of  a  pixel,  and  nu  be  the  desired  boundary 
layer  depth  (in  pixels).  The  steps  are  then  as  follows: 

•  For  i  m  1  to 

-  Create  points  x  =  dj  +  av  —  ipwu ,  a  E  [0,  lj]. 

-  If  p(I(x)\0)  —  p(I(x)\B)  <  0,  then  perturb  the 
points  as  in  (1),  x+(p(I(x)\0)—  p(I(x)\B))gu  — ►  x, 

and  incorporate  them  in  the  data  points  used  to  fit  a 
line  between  vertices  j  and  j  - hi. 


E.  Step  5:  Fit  lines  to  these  perturbed  points  along  each  edge. 

This  is  a  standard  weighted  least  squares  problem  in  the 
plane.  The  weighted  least  squares  fit  to  the  line  uj x  =  bj 
between  vertices  j  and  j  - hi  can  be  found  from  the  points 
along  this  line,  Xk,  by  calculating 


1 

n 


n 

k=  1 


n 

Erp  —  ~T 

wkxkxk  ~nxwxw 

k= l 


where  n  is  the  number  of  perturbed  points  between  vertices 
j  and  j  - hi  and  wk  is  the  weight  on  the  kth  perturbed  point. 
The  unit  normal,  u3 ,  is  then  the  eigenvector  corresponding  to 
the  minimum  eigenvalue  of  M  and 


h3 


AT  _ 
nxwuj 


En 

k=lWk 


(4) 


Since  M  is  only  a  2  x  2  matrix,  the  eigenvalue  problem  can 
be  explicitly  solved  to  speed  the  calculations. 

Points  on  the  object  are  given  a  weight  of  one  (wk  =  1), 
while  any  point  statistical  believed  to  be  on  the  background 
(p(I(x)\0)  —  p(I(x)\B)  <  0)  is  given  much  higher  weight, 
w\).  Our  experiments  indicate  that  at  least  ten  times  higher 
weight  on  background  points  is  useful  for  making  the  lines 
converge  along  the  boundary. 


F.  Step  6:  Create  new  vertices  by  intersecting  non-collinear 
adjacent  lines. 

Let  6  be  the  angle  between  normals  Uj  and  Uj+%.  Then 
Uj+ 1  =  RzUj ,  where 


Rz 


cos  0  —  sin  6 

sin  0  cos  0 


is  a  rotation  matrix  about  the  “z”  axis.  Let 


A?  —  [  A/  A+i  ]  —  [  A/  Rz^j  ] 

The  singular  values  of  Aj  are  by  definition  the  square  root  of 
the  eigenvalues  of  Aj  Aj.  Because  Uj  is  a  unit  vector,  Aj  A  j 
can  be  written  as 


Aj  Aj  =  I  +  uj  Rzitjl1- 


where 


0  1 
1  0 


The  eigenvalues  of  Jx  are  ±1,  thus  the  eigenvalues  of  Aj  A  j 
are  l±uj  RzUj.  Since  u3  is  a  unit  vector  and  Rz  is  a  rotation 
matrix,  uj RzUj  =  cos  0 ,  thus  the  singular  values  of  Aj  are 
y/1  +  cos  0,  s/1  —  cos  9.  The  condition  number  of  Aj  is  the 
ratio  of  the  maximum  and  minimum  singular  values, 


cond(Aj) 


s/l  +  cos  0 
\/l  —  cos  0 


(5) 


For  small  angles  between  adjacent  edges,  Aj  becomes  very 
poorly  conditioned  (cos#  — >  1  =>  cond(Aj)  — >  oo).  As 
the  next  section  will  show,  this  implies  from  (6)  that  vertex 
j  - hi  will  be  poorly  estimated.  To  avoid  this  problem  and 
simultaneously  combine  edges  that  are  nearly  collinear,  cos# 
can  be  found  from  cos  #  =  uj Uj+ \.  If  cos #  <  t,  where  t  is  a 
threshold  less  than  1,  then  a  new  vertex  is  found  by  intersecting 
lines  j  and  j  - hi.  The  new  vertex,  dj+ 1,  is  the  intersection  of 
the  lines  uj x  =  bj  and  uj+1x  =  bj+ \.  Substituting  dj+i  =  x 
into  these  two  equations  and  solving  for  dj+ 1  yields 


dj+i  — 


Uj+i 


T 

=  A~t 

.  kj+i 

3 

.  bj+ 1 

(6) 
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If  the  adjacent  lines  have  nearly  the  same  normal,  but  are 
translated,  then  instead  of  intersecting  the  two  lines,  place  the 
new  vertex  at  the  point  nearest  the  old  vertex  on  the  new  line. 
That  is,  if  cos  6  >  t  and  \bj  —  frj+i|  >  Pu  then  add  the  vertex 
(/  —  Uj+iuJ+1)dj^i  +  bj+\Uj+i.  The  constant  pt  is  a  given 
threshold. 

Because  the  new  list  of  vertices  enforce  collinearity  con¬ 
straints,  the  distance  between  vertices  may  be  much  larger 
than  Id,  and  correspondingly,  the  number  of  vertices  in  the 
new  list  may  be  much  smaller  than  the  number  of  vertices 
found  after  step  3.  This  successive  process  of  expanding,  then 
contracting  the  number  of  vertices  allows  fine  features  to  be 
extracted,  while  simultaneously  automatically  outputting  the 
longest  straight  edges. 

G.  Step  7:  Check  for  convergence,  returning  to  step  2  until 
convergence  takes  place. 

Convergence  can  be  evaluated  by  finding  the  normed  change 
in  vertex  locations.  Once  the  number  of  vertices  has  stabilized, 
they  can  be  stacked  in  a  vector  d  =  [d%  d%  . . .  d^]T - 

The  weighted  norm  of  the  difference  in  this  vector  between 
iterations  is  a  useful  metric  for  determining  convergence.  Until 
it  is  sufficiently  small,  return  to  step  2. 

H.  Establishing  the  Correspondence 

Once  the  straight  line  outline  of  the  object  is  found  using  the 
above  highly  robust  method,  the  correspondences  between  the 
vertices  found  in  the  image  and  those  in  a  CAD  model  can 
be  found.  For  satellites,  a  very  coarse  estimate  of  the  pose 
is  known  from  mission  requirements.  For  instance,  the  solar 
arrays  must  face  the  sun.  This  knowledge  can  be  used  to  find 
which  vertices  should  appear  in  the  image,  and  in  what  order. 

III.  Detection  of  Independently  Moving  but 
Visually  Overlaid  Objects. 

There  exists  a  need  to  detect  an  unwanted  “chaser”  space¬ 
craft  maneuvering  or  orbiting  around  a  satellite.  This  problem 
can  be  broken  into  two  cases: 

1)  The  chaser’s  image  is  distinct  from  that  of  the  satellite. 

2)  The  chaser  and  satellite’s  images  are  blurred  together, 
so  they  are  visually  overlaid. 

In  the  first  case,  it  is  relatively  straightforward  to  detect  the 
presence  of  the  unwanted  chaser  through  a  variety  of  motion 
analysis  techniques.  Consequently,  it  will  not  be  considered 
herein.  The  second  case,  in  contrast,  can  be  difficult:  are  the 
changes  in  the  image  over  time  just  due  to  motion  of  the 
satellite,  or  is  a  chaser  present  and  accounting  for  some  of  the 
blurry  image’s  changes?  The  image  may  be  very  blurry  simply 
because  of  the  distances  involved,  the  available  telescope, 
atmospheric  effects,  or  poor  (perhaps  even  totally  missing) 
illumination.  This  section  will  develop  a  first  level  screening 
procedure  which  can  very  quickly  indicate  the  presence  of 
a  possible  chaser.  Once  this  detection  has  occurred,  further 
resources  can  be  allocated  (more  processing,  better  telescopes, 
human  visual  inspection,  better  illumination  and  atmospheric 
seeing  windows,  etc.). 


The  new  technique  does  not  require  detailed  knowledge 
of  the  satellite’s  shape  or  motion,  but  is  a  simple  method. 
More  complex  algorithms  incorporating  knowledge  of  satellite 
shape  are  derived  in  the  next  section.  From  successive  images, 
this  new,  simple  method  finds  the  dimension  of  the  subspace 
spanned  by  the  motion  in  the  image.  If  the  dimension  is  higher 
than  that  explained  by  the  satellite’s  motion  alone,  then  a 
chaser  may  be  present. 

The  new  method  is  inspired  by  robotic  motion  segmentation 
algorithms  [18].  Suppose  the  satellite  has  n  points  that  are 
seen  on  all  of  the  F  image  frames.  Let  the  jth  satellite 
point  be  denoted  as  pj  (j  =  1 . . .  n).  Let  g  be  the  4  x  4 
homogeneous  coordinate  transformation  from  the  satellite’s 
coordinate  system,  S,  to  the  camera’s  coordinate  system,  C. 
Then  [5] 


'  cpj  ' 

'  sPj ' 

1 

—  9 

i 

where  the  prescripts  C  and  S  denote  that  the  point  is  measured 
in  the  C  and  S  coordinate  systems,  respectively.  Denoting 


equation  (7)  can  be  written  more  compactly  as 

CPj  =  9S  Pj  (8) 

Key  fact:  At  a  given  time,  the  same  g  is  used  to  transform 
all  points  on  the  satellite  into  camera  coordinates. 

Let  the  camera  have  a  focal  length  of  /,  and  let  the  satellite 
be  at  a  distance  z.  If  the  camera’s  axis  is  the  z  axis,  then  the 
image  plane  position  of  pj  in  the  ith  frame  is 

da  =  ~z  [  h  o  ]  cpj 

Note  that  this  equation  assumes  that  the  distance  from  the 
camera  to  the  satellite  is  much  larger  than  the  size  of  the 
satellite,  so  minor  differences  in  the  distance  from  the  camera 
to  a  point  on  a  satellite  caused  by  that  point’s  location  on 
the  satellite  are  irrelevant  (affine  camera  model).  The  affine 
camera  model  is  highly  accurate  when  imaging  objects  from 
Earth  to  space,  as  the  distance  is  hundreds  of  kilometers  even 
for  low  Earth  orbits.  From  (8), 

dij  =  {  [  h  0  ]  gfpj  (9) 

where  the  i  subscript  on  g  denotes  the  homogeneous  coordi¬ 
nate  transformation  when  the  ith  frame  is  taken.  Denote  the 
2x4  transformation  matrix  as 

Mi  =  -  [  I2  0  ]  gi 
z  L  J 

Equation  (9)  then  becomes 

dij  =  Mi  pj 

Let  the  n  points  be  tracked  from  one  frame  to  the  next  over 
i  =  1 ...  F  frames.  Then  the  ith  row  of  the  data  matrix,  D  = 
{d^}  is 

Mi  [5P!  sp2  ■■■spn  } 
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Denote  the  4  x  n  post-multiplying  matrix  as  SP , 


SP=  [SP1  SP2  SPU  ] 


Finally, 


D  = 


Mi 

M2 


Mf 


2Fx4 


4xn 


(10) 


Equation  (10)  shows  that  the  data  matrix,  D  =  {dij}  can 
be  decomposed  into  two  matrices  multiplied  together.  Since 
image  data  is  two  dimensional,  the  first  matrix  has  dimension 
2 F  x  4,  while  SP  is  4  x  n.  As  the  rank  of  a  matrix  is  less  than  or 
equal  to  the  smallest  dimension  of  its  factors,  the  rank  of  D  is 
less  than  or  equal  to  four.  This  is  a  surprising  result,  as  D  may 
be  a  very  large  matrix  if  many  points  are  tracked  across  many 
frames.  If  only  rotations  are  present,  then  a  similar  analysis 
can  be  used  to  show  that  the  rank  of  D  is  less  than  or  equal 
to  three. 

When  the  real  data  is  collected,  some  of  the  points  may  not 
in  actuality  be  located  on  the  satellite-they  could,  for  instance, 
lie  on  a  chaser  vehicle.  In  this  case,  the  above  arguments  no 
longer  hold,  and  the  dimension  of  D  can  exceed  four.  This 
implies  that  the  dimension  of  the  data  matrix  can  be  used  as  a 
detection  criteria-if  it  exceeds  four,  then  a  chaser  may  possibly 
be  present,  so  further  investigations  are  merited. 

Note  that  this  detection  algorithm  is  very  simple  to  imple¬ 
ment,  as  it  requires  no  knowledge  of  the  satellite.  Moreover, 
point  correspondences  are  only  required  between  frames.  This 
is  comparatively  simple,  since  image  points  typically  don’t 
move  much  over  a  single  sample  time,  thus  finding  two  corre¬ 
sponding  points  is  simplified.  To  implement  the  algorithm,  the 
measurements  must  simply  be  collected  in  a  single  matrix,  D , 
where  each  column  contains  the  image  positions  of  a  single 
satellite  point  as  it  appears  in  different  images.  The  rank  of  D 
is  calculated  to  detect  multiple  moving  bodies. 


IV.  Detection  of  Extraneous  or  Missing  Objects 

The  simple  technique  in  Section  III  is  very  useful  for 
many  scenarios,  but  it  can  be  improved  by  including  more 
knowledge.  This  section  assumes  that  a  solid  model  of  the 
satellite’s  shape  is  available,  and  uses  that  knowledge  to 
compare  the  obtained  image  with  predictions  of  what  the 
image  should  look  like.  Gross  differences  indicate  an  anomaly. 

This  section  briefly  reviews  pose  estimation  results  found 
in  [12].  The  proofs  and  further  details  (including  methods  for 
evaluating  the  estimation  error  and  viewpoint)  are  in  [12]. 
From  a  noisy,  diffraction  limited  image  of  a  known  satellite 
in  an  unknown  pose,  the  pose  will  be  estimated.  The  CAD 
model  can  then  be  placed  in  the  same  pose  and  projected  onto 
two  dimensions  to  create  a  synthetic  prediction  of  what  the 
satellite’s  image  should  look  like.  Subtraction  of  the  actual 
image  from  the  synthetic  image  produces  a  signal  useful  for 
automatically  checking  any  satellite  anomalies. 


First,  the  estimation  of  orientation  alone  will  be  treated. 
Define  the  hat  function  (•)  as  the  cross  product  matrix,  i.e. 


0 

-Z3 

Z2 

£3 

0 

-zi 

-Z2 

zi 

0 

Theorem  1:  Given  Zi,  Vi  E  M3  and  Wi  E  R+,  i  —  1 . . .  n, 
the  minimum  of 


n 

J  =  Wj\\Rzj  —  avi  ||2 


over  R  E  SO(3),  a  E  M+  is  found  by  first  calculating 

n 

D  =  'P,wiK(yi,Zi) 
i=  1 

where  K  is  the  bilinear,  symmetric  matrix  function 


K{v,z) 


VZ1 


+  ziF  —  2  vFzI  zv 
(zv)T  0 


(11) 


(12) 


(13) 


The  eigenvector,  e,  of  D  corresponding  to  the  maximum 
eigenvalue  is  the  unit  quaternion  representation  of  the  optimal 
R.  The  optimal  scaling  is 


=  ELiMM 

"  E"=i  ^11^1 12 


(14) 

□ 


Orientation  estimation  from  2-d  data  does  not  enjoy  the 
invariance  to  scaling  found  in  Theorem  1,  so  the  case  where 
scaling  is  known  will  be  solved,  then  scaling  will  be  estimated 
separately. 

Proposition  2:  Given  zi  E  M3,  di  E  M2,  and  Wi  E  R+, 
i  =  1 . . .  n,  P  =  [I 2  0]  the  minimum  of 

n 

J2d  =  ^WiWPRzj  -  dj ||2  (15) 

i=i 


over  R  E  SO(3),  is  found  by  first  calculating 


Vi  = 


±vW-||  di 


,  i  =  1 . . .  n 


(16) 


Next,  the  2n  3-d  orientation  problems  corresponding  to  all 
possible  Vi ,  i  =  1 . . .  n  d=  sign  choices  are  solved  using 
Theorem  1: 

n 


min  >  wA\Rz 
Reso(  3)^ 

i=i 


i 


From  these  2n  solutions,  the  optimum  is  that  minimizing  J^d- 


□ 


Lemma  3  will  now  show  how  these  results  can  be  brought 
together  to  estimate  orientation  from  2-d  measurements  with 
unknown  scaling. 

Lemma  3:  Given  Zi  E  M3,  di  E  M2,  and  Wi  E  M+,  i  = 

1 . . .  n,  P  =  [I2  0]  the  minimizer  of 

n 

J  =  ^Wi\\PRzi~  adi\\2  (17) 

i=  1 
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over  R  G  SO(3),  a  G  M+,  is  approximated  by  first  estimating 
ft  by  ae  =  - ~t  =r-  77—  - .  R  can  then  be  estimated 

J  e  maxl=i...n  \\di\\/\\ziW 

by  replacing  all  measurement  vectors  di  with  scaled  versions 
di  1  *  aedi,  i  =  1 . . .  n,  and  applying  Prop.  2. 

n 


Finally,  to  estimate  both  orientation  and  position,  use  Prop. 
4. 

Proposition  4:  Given  data  points  zt  G  M3  and  correspond¬ 
ing  sensed  points  qi  G  M2,  along  with  weightings  Wi  G  M+, 
i  =  1 . . .  n,  P  =  [1 2  0] ,  the  minimizer  of 


n 

J  =  YJWi\\P[Rzi+p\-aqi\\2  (18) 

i= 1 


over  g  =  (R,p)  G  SE(3)  ( R  G  SO(3),  p  G  M3),  and  a  G  M+ 
can  be  estimated  by  forming  the  free  vectors  %  =  zi  —  Zj , 
=  Qi  ~  Qj>  h  j  =  1  •  •  •  n.  and  a  can  then  be  estimated 
using  Lemma  3  by  replacing  Zi  by  and  d{  by  d^ .  The  first 
two  components  of  p  are  estimated  by 


Px 

_  Py  _ 


=  Pa  =  PP  = 


a  Y!i= 1  wi4i  ~  Ra  Y!i= 1  wi*i 

En 

i=lWi 


(19) 


where  Ra  denotes  the  2  x  3  matrix  consisting  of  the  first  two 
rows  of  R.  When  the  measurements  arise  from  an  optical  sys¬ 
tem  with  focal  length  /,  and  pz  denotes  the  third  component 
of  p,  and  if  the  imaging  is  affine  (||^||  «])2  for  all  points 
i  =  1 . . .  n),  then 

Pz  ~  af  (20) 


□ 


V.  Simulation  Results 

Figures  1  and  2  illustrate  the  silhouettes  obtained  using  a 
standard  gradient  (Sobel)  operator  for  both  the  sun  illuminated 
and  un-illuminated  (only  back  lit)  cases.  While  the  Sobel 
operator  performs  well  when  the  satellite  is  illuminated,  it 
completely  fails  without  illumination. 

The  new  SSLS  method,  on  the  other  hand,  is  able  to  find 
the  silhouette  (Fig.  3).  Table  4  documents  the  parameters  used 
in  this  simulation. 

Figure  5  illustrates  a  sequence  of  diffracted  images  con¬ 
taining  both  the  satellite  and  a  chaser  in  a  nearby  relative 
orbit.  The  +  signs  indicate  points  that  are  tracked  between 
each  image,  while  the  circled  +  signs  indicate  those  points 
due  to  the  chase  vehicle.  The  satellite’s  translational  motion  is 
tracked,  so  only  rotations  are  present  in  the  image.  A  data 
matrix,  D ,  is  formed  where  each  column  lists  the  image 
plane  positions  of  a  single  point.  Since  only  rotations  are 
present,  a  single  rigid  body  ideally  has  rank(D)  of  at  most 
three.  When  all  eight  columns  are  included,  the  singular 
values  are:  (347.28,75.33,31.48,3.35,0.69,0.51,0.23,0.14). 
If  the  two  columns  corresponding  to  points  on  the 
chaser  are  excluded,  then  the  singular  values  become: 
(307.55,73.49,1.32,0.63,0.48,0.19).  In  the  later  case,  the 
rank  should  be  at  most  three,  since  all  the  data  points  stem 
from  a  single  rigid  body.  However,  the  rank  is  greater  than 
three  in  this  realistic  simulation  due  to  noise  and  imperfect 


Illuminated  Satellite 
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Fig.  1.  The  illuminated  satellite  (top),  and  its  silhouette  found  using  gradient 
masks  (bottom). 


correspondences.  Nevertheless,  the  third  and  higher  singular 
values  in  the  multiple  body  case  are  significantly  larger  than 
those  in  the  later,  single  body  case-a  discernable  difference 
that  can  be  wisely  used  as  an  indication  of  multiple  body 
presence. 

From  a  noisy,  diffraction  limited  image  of  a  known  satellite 
in  an  unknown  pose,  the  pose  will  be  estimated.  The  CAD 
model  can  then  be  placed  in  the  same  pose  and  projected  onto 
two  dimensions  to  create  a  synthetic  prediction  of  what  the 
satellite’s  image  should  look  like.  Subtraction  of  the  actual 
image  from  the  synthetic  image  produces  a  signal  useful  for 
automatically  checking  any  satellite  anomalies. 

Fig.  6  depicts  a  hypothetical  satellite.  The  far  corners  of  the 
solar  arrays  are  chosen  as  the  data  points  (zi,  i  =  1, ...  ,4, 
depicted  as  0’s),  as  they  produce  an  easily  identifiable  signa¬ 
ture.  As  the  satellite  moves,  several  images  are  captured  (Fig. 
7).  The  pose  will  be  estimated  and  the  satellite  inspected  for 
images  1  and  4  (1  is  the  leftmost). 

Because  satellites  orbit  at  significant  distances  from  the 
Earth,  the  actual  image  is  degraded  appreciably  by  diffraction 
effects  (see  Fig.  8).  Using  a  binary  thresholding  technique,  it 
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Un-illuminated  Satellite  with  Parasite 
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Fig.  2.  The  un-illuminated  satellite  image  (top),  and  its  gradient  (bottom). 
The  silhouette  is  not  found. 

Silhouette  of  Un-llluminated  Satellite:  Straight  Line  Snakes 
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Silhouette  of  Un-llluminated  Satellite:  Sobel  Operators 
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Fig.  3.  The  silhouette  is  found  using  the  SSLS  method. 
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Fig.  4.  Parameters  used  in  the  SSLS  method. 


Fig.  5.  A  sequence  of  diffracted  images,  and  the  tracked  correspondence 
points. 
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Fig.  6.  CAD  model  of  the  imaged  satellite,  with  O’s  denoting  the  point 
correspondences. 


Sequence  of  Satellite  Images. 
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Fig.  7.  As  the  satellite  sweeps  across  the  sky,  several  images  are  captured. 

is  enhanced,  and  the  solar  array  corners  in  the  image  ($)  are 
found  (Fig.  9). 

Despite  the  effects  of  diffraction,  thresholding,  and  pixel 
round-off,  the  pose  estimate  is  quite  accurate,  producing  an 
angular  estimation  error  of  only  eg  =  5.3°.  A  thresholded  ver¬ 
sion  of  the  predicted  image  is  subtracted  from  the  thresholded 
measured  image  (Fig.  9)  to  produce  the  inspection  image  (Fig. 
10).  Although  it  suffices  for  our  present  purposes,  an  improved 
match  can  certainly  be  obtained  from  this  subtracted  image  by 
employing,  for  example,  the  methods  in  [7]  to  fine  tune  the 
estimate. 

Automatically  determining  the  point  correspondences  (fo  can 
be  difficult.  These  simulations  first  find  a  point  inside  the 
satellite  by  a  very  coarse  template  matching  strategy.  Then, 
a  statistical  snake  [14]  grown  from  that  initial  point  robustly 
identifies  the  perimeter.  Finally,  vertices  along  that  perimeter 
are  extracted  (see  [13]  for  further  details).  Figures  11  and  12 
illustrate  the  results  on  image  4.  In  this  case,  eg  =  25°. 

VI.  Conclusions 

This  paper  develops  three  new,  interconnected  techniques 
useful  for  the  autonomous  distant  visual  surveillance  of  satel¬ 
lites.  First,  silhouetting  of  man  made,  erratically  illuminated 


Diffracted  Image 
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Fig.  8.  The  first  raw,  diffracted  image. 


Satellite  Image  and  its  Estimated  Vertices  (o) 
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Fig.  9.  The  first  image  with  the  sensed  points,  qi,  denoted  by  (0). 


Matched  Image  minus  Measured — gray  means  match  is  good 
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Fig.  10.  Combining  pose  estimates  with  the  CAD  model,  a  predicted  image 
is  calculated  and  subtracted  from  the  measured  image. 
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Satellite  Image  and  its  Estimated  Vertices  (o) 
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Fig.  11.  The  fourth  image,  with  its  points  cfo  denoted  by  (0).  The  outside 
perimeter  is  identified  by  statistical  snakes. 


Matched  Image  minus  Measured — gray  means  match  is  good 


Fig.  12.  Combining  pose  estimates  with  the  CAD  model,  a  predicted  image 
is  calculated  and  subtracted  from  the  fourth  measured  image. 


satellites  is  performed.  Illumination  cases  include  full  sun  from 
an  arbitrary  (often  awkward)  viewing  angle  and  un-illuminated 
(back-lit  by  the  star  field).  New  Statistical  Straight  Line  Snakes 
(SSLS)  prove  efficient  in  finding  the  silhouette,  even  in  the 
un-illuminated  case.  In  contrast,  gradient  techniques  based  on 
the  Sobel  operator  fail  utterly.  The  silhouette  is  composed 
of  straight  line  segments,  which  are  easy  to  calculate,  fit 
the  straight  lines  inherent  in  man  made  objects,  and  lend 
themselves  to  further  processing  (pose  estimation,  template 
matching,  etc.)  A  current  silhouette  can  be  used  as  an  initial 
guess  to  speed  silhouette  tracking  between  successive  images. 
Once  the  silhouette  has  been  used  to  find  correspondence 
points,  a  second  method  for  detecting  a  moving,  nearby 
chaser  vehicle  is  derived.  The  hard  case  is  treated  in  which 
the  chaser  and  satellite  are  so  nearby  that  their  images  are 
blurred  together.  The  algorithm  finds  the  dimension  of  motion 
generated  by  the  sequence  of  images.  If  the  dimension  is 


higher  than  that  explained  by  a  single  rigid  body,  then  this 
indicates  a  possible  chaser.  No  knowledge  of  the  satellite  or 
chaser’s  shape  or  motion  is  required.  Simulations  indicate 
discemable  differences  in  the  singular  values  obtained  with 
one  versus  two  bodies  present. 

Independent  relative  motion  between  the  satellite  and  chaser 
is  required-if  the  chaser  is  immobile  with  respect  to  the 
satellite,  then  a  third  technique  must  be  used.  This  third 
method  incorporates  the  satellite’s  solid  model  to  estimate 
its  pose  from  a  noisy,  diffraction  limited  image.  The  pose 
is  then  combined  with  the  solid  and  optical  model  to  create 
synthetic  expected  images.  Inspection  is  performed  by  com¬ 
paring  these  with  the  actual  images.  The  new  pose  algorithm 
first  estimates  depth  by  a  least  upper  bound  technique.  A 
fast  method  is  derived  of  optimally  estimating  the  rotation 
matrix  by  a  sequence  of  analytical  solutions  (rather  than  a 
nonlinear  numerical  optimization!).  Simulations  illustrate  a 
satellite  inspection  from  Earth  to  low  Earth  orbit. 
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